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Abstract. Simulation results of dense granular gases with particles of different size 
are compared with theoretical predictions concerning the pair-correlation functions, 
the coUison rate, the energy dissipation, and the mixture pressure. The effective 
particle-particle correlation function, which enters the equation of state in the same 
way as the correlation function of monodisperse granular gases, depends only on 
the total volume fraction and on the dimensionless width A of the size-distribution 
function. The global equation of state is proposed, which unifies both the dilute and 
the dense regime. 

The knowledge about a global equation of state is applied to steady-state situ- 
ations of granular gases in the gravitational field, where averages over many snap- 
shots are possible. The numerical results on the density profile agree perfectly with 
the predictions based on the global equation of state, for monodisperse situations. 
In the bi- or polydisperse cases, segregation occurs with the heavy particles at the 
bottom. 



1 Introduction 

The hard-sphere (HS) gas is a traditional, simple, tractable model for various 
phenomena and systems like e.g. disorder-order transitions, the glass transi- 
tion, or simple gases and liquids 0-|3|. A theory that describes the behavior 
of rigid particles is the kinetic theory |l],^, where particles are assumed to 
be rigid and collisions take place in zero time (they are instantaneous), ex- 
actly like in the hard-sphere model. When dissipation is added to the HS 
model, one has the most simple version of a granular gas, i.e. the inelastic 
hard sphere (IHS) model. Granular media represent the more general class of 
dissipative, non-equilibrium, multi-particle systems Attempts to describe 
granular media by means of kinetic theory are usually restricted to certain 
limits like constant or small density ||] or weak dissipation [0,^. Also in the 
case of granular media, one has to apply higher order corrections to success- 
fully describe the system under more general conditions Already clas- 
sical numerical studies showed that the equation of state can be expressed as 
some power series of the density in the low density regime [^,|^, whereas, 
in the high density case, the free volume theory leads to useful results |jl^. In 
the general situation, the granular system consists of particles with different 
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sizes, a situation which is rarely addressed theoretically p^-p^]. However, the 
treatment of bi- and polydisperse mixtures is easily performed by means of 
numerical simulations JlTj-p^. 

In this study, theories and simulations for situations with particles of 
equal and different sizes are compared. In section || the model system is 
introduced and in ^ we review theoretical results and compare them with 
numerical results concerning correlations, collision rates, energy dissipation 
and pressure. Based on the numerical data, a global equation of state is 
proposed. This global equation of state is valid for arbitrary densities and 
mixtures with particles of different size, and it is used to explain the density 
profile in a dense system in the gravitational field in section ^. The results 
are summarized and discussed in section |^. 

2 Model system 

For the numerical modeling of the system, periodic, two-dimensional (2D) 
systems of volume V = L^Ly are used, with horizontal and vertical size 
and Ly, respectively. N particles are located at positions with velocities 
Vi and masses m.^. From any simulation, one can extract the kinetic energy 
E = ^ '^i'^'ii dependent on time via the particle velocity Vi. In 2D, the 
"granular temperature" is defined as T = E/N . 

2.1 Polydispersity 

The particles in the system have the radii a, randomly drawn from size dis- 
tribution functions wia) as summarized in table |l| where the step- function 
9[x\ = 1 for a; > 1 and 9[x\ = for a; < 1 is implied. 



(i) 


monodisperse 


w{a) — S{a — ao) 


(ii) 


bidisperse 


w{a) — niS{a — ai) + n2(5(a — 0,2) 


(iii) 


polydisperse 


2u,nan^[" (1 wo)ao]e[{l + wo)ao a] 



Table 1. Size distribution functions used in this study. 



The parameter og is the mean particle radius (a) in cases (i) and (iii) . In the 
bidisperse situation (ii), one has Uq = (a) = niOi-l- 7^202 = (ni-|-(l— 7ii)/i?)ai, 
with the fraction rti = Ni/{Ni + N2) of particles with size ai in a system 
with A^i + A^2 particles in total and N2 particles with radius 02 . Thus, 
besides ni, only the size ratio R = axja-i is needed to classify a bidisperse 
size distribution. The total volume fraction v — v\ ^ V2 Ss the last relevant 
system parameter, since the partial volume fractions z^i,2 — N\,2T^a\ 2/^ — 
nx.iva^ 2/ ('^^) can be expressed in terms of n\ and R: Using the dimensionless 
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moments 

A, = m + (1 - = ^ , (1) 

one has i^i — niv/A2 and V2 = {1 — ni)v/{E?A2). Since needed later on, the 
expectation values for the moments of a and their combination, the dimen- 
sionless width-correction A = (a)^/ (a^), are summarized in table ^ in terms 
of ai, and R for the bidisperse situations and in terms of ao and wq in the 
polydisperse cases. Different values of v are realized by shrinking or growing 
either the system or the particles. 







(a) 




(«)7K> 


(i) 


monodisperse 


ao 


ao 


1 


(ii) 


bidisperse 




A2ai 


Ai/A2 


(iii) 


polydisperse 


ao 


[l + wl/Z] at, 


3 / (3 + ™c1) 



Table 2. Moments (a), {a?) and A = (a)^ / {a?) of the size distribution functions. 



2.2 Particle Interactions 

The particles are assumed to be perfectly rigid and to follow an undisturbed 
motion until a collision occurs as described below. Due to the rigidity, col- 
lisions occur instantaneously, so that an event driven simulation method 
[^|T| can be used. Note that no multi-particle contacts can occur in this 
model. For a review on possible, more physical extensions of this model see 
Ref. H]. 

A change in velocity - and thus a change in energy - can occur only 
at a collision. The standard interaction model for instantaneous collisions of 
particles with radii a.i, mass rrii = (4/3)7r/3a?, and material density p is used 
in the following. (Using the mass of a sphere is an arbitrary choice, however, 
using disks would not influence most of the results discussed below.) This 
model was introduced and also discussed for the more general case of rough 
particle surfaces in Refs. [[7| |2^-|25|] . The post-collisional velocities v' of two 
collision partners in their center of mass reference frame are given, in terms 
of the pre-coUisional velocities by 

, (1 + "^12 
V'l 2 = ^'1,2 T ^ «n , 2 

with Vn = \{vi — V2) ■ n] n, the normal component of Vi — V2 parallel to n, 
the unit vector pointing along the line connecting the centers of the colliding 
particles, and the reduced mass mi2 = mi77i2/(mi -I- m2). If two particles 
collide, their velocities are changed according to Eq. (||) and any r(v„), de- 
pendent on ti„ = \vn\, can be used. For a pair of particles, the change of the 
translational energy at a collision is AE = — mi2(l — r^)w^/2. 
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3 Simulation and theory 

In the following, we compare simulations with different polydispersity, i.e. dif- 
ferent size distribution functions w{a), as summarized in table ^. 







w{a) parameters 


particles 


A 


A 


moriodisperse 


Wo = 0, 


iV = 1628 


1 


B 


monodisperse 


ni = 1, _R = 1 


iV = 576 


1 


C 


bidisperse 


m = 0.517, i? = 3/4 


iV = 576 


0.9798 


Dl 
D2 


bidisperse 
bidisperse 


m = 0.781, i?= 1/2 
ni = 0.799, i?= 1/2 


iV = 576 
iV = 6561 


0.8968 
0.8998 


E 


polydisperse 


wo = 0.25 


iV = 1425 


0.9796 


Fl 

F2 


polydisperse 
polydisperse 


Wo = 0.5 
Wo — 0.5 


iV = 1425 
iV = 1521 


0.9231 
0.9231 



Table 3. Simulation parameters for the simulations discussed below. Note that sets 
C and E have different w{a) but almost identical values of A. 



3.1 Particle correlations 

In monodisperse systems, the particle-particle pair correlation function at 
contact, 

1 - 7u/16 

92a{v) - _ , (3) 

depends on the volume fraction only [0,^,|^,||,|2^. The particle-particle cor- 
relation function is obtained from the simulations by averaging over M snap- 
shots with N particles, normalized to the value g{r ^ 2a) = 1 for long 
distances in large systems, so that 

M N 

= M E W^^V E E ^[-^. - + ^r- n,] , (4) 

m=l ^ ' ^ i=l j=l 

with rij — jr — rj I, and where the two functions select all particle pairs (i, j) 
with distance between r and r + Ar. The weight N{N — l)/2 accounts for all 
pairs summed over, and the term Vr = 7r(2r -I- Ar)Ar is the volume (area) of 
a ring with inner radius r and width Ar. In Fig. |l], simulation results from 
set B (see table |3|) with different v are presented. Typical values used for the 
averages are e.g. M = 50, and Ar = a/10. Besides fluctuations, the values at 
contact nicely agree with the theoretical predictions from Eq. (||) , as indicated 
by the arrows - as long as the system is disordered (left panel in Fig. |^). In 
a more ordered system (right panel in Fig. |^), g2a{v) is not a good estimate. 
Instead, one obtains a long range order with peaks at r/2a = l,\/3, 2, ... , 
indicating the triangular lattice structure of the assembly. 
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Fig. 1. Particle-particle correlation function g{r) plotted against the normalized 
center-center distance r/2a. (Left) Disordered systems - the arrows indicate the 
values at contact from Eq. (Right) Ordered systems with different axis scaling. 




Fig. 2. Particle-particle correlation function g{r) plotted against the center-center 
distance normalized by the radius ai of the smaller particles, at a volume fraction 
of v = 0.60. The arrows indicate the values gn, gi2, and g22 at contact. (Left) 
Bidisperse simulation from set C with R = 3/4, and (Right) bidisperse simulation 
from set Dl with R=l/2. 



For bidisperse situations, the pair correlation functions for equal species 
gii and 522 are obtained by replacing N in the first particle-sum in Eq. 
by iVi and N2, respectively. For the correlation function 1712, it is necessary 
to perform the first particle-sum in Eq. (|^) from i — 1 to Ni, the second sum 
from j = 1 to N2, and to replace the weight N{N — l)/2 by N1N2 (in order 
to account for all pairs of different kind). In Fig. ^ simulation results from 
sets C and Dl for i/ = 0.6 are compared to the analytical expressions Eqs. 




Fig. 3. Particle-particle correlation function g{r) plotted against the normalized 
center-center distance r/2ao. (Left) Low density systems and (Right) high density 
systems with the same axis scaling. 



(90), (91), and (92) from Ref. [Q, here expressed in terms of Ai_2, -R, and v. 



.911 = n — ' (5) 



16 A2 



9 Ax 

522 = 7^ , and (6) 



1-1^)2 

1 

(T 



512 - T. rT2 ■ (7) 



Note that all are identical to gia(y) in the monodisperse case with i?, = 1 
and A\ = = 1. The parameters used for averaging were M = 50 and 
Ar = ai/118 (Left) and Ar = ai/64 (Right). A finer binning leads to 
stronger fluctuations, a rougher binning does not resolve the values at con- 
tact, however, within the statistical error, the agreement between theoretical 
predictions and numerical results is reasonable. 

Finally, in Fig. ^, particle correlation functions from polydisperse simula- 
tions (set Fl) are presented. Due to the broad and continuous size distribu- 
tion function, g(r) is a smooth function with much less variety in magnitude 
than in the mono- and polydisperse situations discussed above. It more re- 
sembles the distribution function of a gas or liquid with a smooth interaction 
potential [||. 



3.2 Collision rates and energy dissipation 

In order to estimate the rate of change of energy in the system T = t~^^AT, 
the collision frequency is needed. Rather than going into details con- 



Poly disperse granular gases 7 



cerning the calculation of t^^, we will simply use the Enskog collision rate 
[|[^,|2^ for identical particles, 

= -^^9ia(v)^jT jm , (8) 
and, equivalently, the inter-species collision rates 

where all rates give the number of collisions of a particle per unit time, with 
ciij = tti + ttj. The temperature is here assumed to be independent of the 
particle species for the sake of simplicity. This is approximately true in the 
systems examined below, provided they stay rather homogeneous, but it is 
not true in general, since the cooling rates depend on the species. 

In Eq. d), the term ^r/(my) is proportional to the mean relative ve- 
locity vlf of a pair (i, j), so that the remainder tij^/vlf can be seen as a 
measure for the inverse interspecies mean free path Ay . The mean collision 
rate in the system is 

^n7ix = Yl " ^^VTrV^^/wi X] ^i^jdijCij , (10) 

with cii ^ 1, ci2 = (1 + R)/{2R)y/{l + R^)/2, and C22 — VR- Note that 
ci2 and C22 depend on mass and density of the different species. The mean 
collision rate was tested for the monodisperse and bidisperse situations and 
showed the same quality of agreement as the pressure, which will be discussed 
in the next subsection. Therefore, we do not present data of the collision 
rates here, but perform a detailed numerical study of the mixture pressure 
below. However, we should remark that the interspecies collision rates are 
of the same order of magnitude, even if the species fluctuation velocities 
Vi = \JT jmi strongly differ due to the differences in mass, also when the 
temperature T = Ti — Tj — E/N is not species dependent. This means 
that the mean free distance between collisions compensates the speed; the 
distance traveled between collisions is proportional to the species velocity. In 
Fig. H several snapshots from simulation set Dl are presented. The grey-scale 
indicates the collision rates, dark particles collide more frequently than light 
particles, and the collision rate increases with the density, but as discussed 
above, the collision rates of the two species are comparable. 

Knowing both the collision rates and the energy loss per collision 

AT,, = -m.,^-^{vlfr - -^-^T , (11) 
it is straightforward to compute the decay of energy as a function of time 

^ = E ^^t^^E,, = —T^-lit)Tit) , (12) 
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V = 0.30 



V = 0.60 




V = 0.70 



ly = 0.85 




Fig. 4. Snapshots from the bisdisperse simulations (set Dl) for different volume 
fractions u. The grey-scale denotes the collision rate of the corresponding particle, 
the darkest particles had a collision rate of 1500 s~^, the lightest particles had a rate 
of less than 275s~^. (These numbers have to be multiplied with 2/5 for v = 0.30 
and with 4 for v = 0.85, in order to allow for a comparison of all pictures.) 



where both T(t) and T[^-^it) oc \/T(t) depend on time. The differential equa- 
tion is easily solved and one gets the scaled temperature 

^ = (i + ^^^,;.i(o)t) ' , (13) 

identical to the solution of the homogeneous cooling state of monodisperse 
disks m, where t^^(O) is replaced by T^^iliO) from Eq. 

In Fig. ^, simulations from set D2 are presented for different values of 



< 1. The agreement between simulations and Eq. (13) is perfect for short 
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times. With decreasing r, i.e. increasing dissipation, the deviations from the 
theory occur earlier due to the break-down of the homogeneity and the related 
simplifying assumptions of molecular chaos and Gaussian velocity distribu- 
tions |,|2|]. 




r=0.999 
r=0.998 
r=0.995 
r=0.990 
r=0.980 
r=0.950 
r=0.900 
r=0.850 



Fig. 5. Dimensionless temperature T /Tq plotted against rescaled time r = ^,7ix(0)^ 
for different r, with double logarithmic axis. The symbols are simulation results 
from set D2, the solid lines correspond to Eq. (p^). In the inset the quality factor 
qT = rsim/Ttheory IS plotted agaiust the time t. 



3.3 Stress and the equation of state 

The stress tensor, defined for a test- volume V , has two contributions, one from 
the collisions and the other from the translational motion of the particles. 
Using a and b as indices for the cartesian coordinates one has the components 
of the stress tensor (where the sign is convention) 



V 



1 



(14) 



with i'j, the components of the vector from the center of mass of the two col- 
liding particles j to their contact points at collision n, where the momentum 
Apj is exchanged. The sum in the left term runs over all particles i, the first 
sum in the right term runs over all collisions n occuring in the time-interval 
At, and the second sum in the right term concerns the collision partners of 
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collision n - in any case the corresponding particles must be within the aver- 
aging volume V pl|,[2^-[2^. Note that the results may depend on the choice 
of At and V |2^, however, a discussion of different averaging procedures and 
parameters is far from the scope of this study. 

The mean pressure p — (cti +cr2)/2 , with the eigenvalues cri and <T2 of the 
stress tensor, can be obtained from simulations with rigid, elastic particles 
(r = 1) and different volume fractions v |^,|2^. The dimensionless reduced 
pressure from simulations agrees perfectly with the theoretical prediction 

PQ^PV/E-l = 2vg2a{v) , (15) 

with the total energy E — {l/'l)E^^mivf. In Fig. || the simulation results 
for the dimensionless pressure P — pV/ E — 1 are compared to the kinetic 



theory result Pq = 2z^52a('^) 21 



100 



10 : 



0.1 



Wo=0.00 
Wo=0.25 
Wo=0.50 
/?=0.50 
R=Q.15 
P, 




0.2 0.4 0.6 

J I 



0.2 



0.4 0.6 
V 



0.8 



Fig. 6. Dimensionless pressure P = pV/E— 1 plotted against the volume fraction v 
for different particle size-distribution functions. The bidisperse situation is identified 
by the size ratio R. The solid line corresponds to Pq. In the inset the quality factor 
P/Po is plotted for the same data. The wiggle of the wq — data is discussed in 
more detail below. 



When plotting P against the volume fraction u with a logarithmic vertical 
axis, the results for the different simulations can not be distinguished for v < 
0.7. In the monodisperse system, we obtain crystallization around ly ~ 0.7, 
and the data clearly deviate from Pq, i.e. the pressure is strongly reduced 
due to crystallization and, thus, enhanced free volume. The monodisperse 
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data diverge at the maximum packing fraction t'™™" — 7''/(2\/3) in 2D. 
Note that one has to choose the system size such that a triangular lattice 
fits perfectly in the system, i.e. Ly/L^ = V3h/2w with integer h and w 
— otherwise the maximum volume fraction is smaller. All other simulations 
are close to Pq up to w 0.8, where they begin to diverge; the bidisperse 
data with R = 1/2, for example, diverge at v^^^ ~ 0.858. The maximum 
packing fraction is smaller for the polydisperse size distributions used here 
and the crystallization, i.e. the pressure drop, does not occur for polydisperse 
packings with wq ~ 0.15 (the data which lead to this approximate result 
are not shown). Since the logarithmic axis hides small deviations, we plot 
also the quality factor P/Pq of the data in the inset. In this representation, 
values of unity mean perfect agreement, while smaller values correspond to 
an overestimation of the data by a factor of Pq/P when Pq would be used 
instead of the simulation results P. The deviations increase with increasing 
width of w{a) and with increasing volume fraction. Note that there exists a 
deviation already for small v. 

A more elaborate calculation in the style of Jenkins and Mancini, see 
Eq. (60) in |l^, leads to the partial translational pressures p* = UiE/V for 
species i and to the coUisional pressures p'ij = nNiNjgijafj{l + rij)T / {AV^) 
with the particle correlation functions from Eqs. (||)-(0) evaluated at contact, 
and aij = + aj. In the simulations from Fig. |^, the inter-species restitution 
coefficients are equal and elasticity is assumed, r = rn = ri2 = r22 = 1- Note 
that the species temperatures are equal, so that the corresponding correction 
term can be dropped. Thus, the global pressure in the mixture is 



1 + (1 + r) ^ ^ {giia\^n\ + 25i2af2niri2 +32202 

E 



A2afi 



= -[l + {l + r)vgj,{v)] . (16) 
Assuming a monodisperse system as a test case, i.e. inserting R = Ai ^ 



A2 = I, into Eq. (|16D, leads to the monodisperse solution p™V/E — 1 = Pq, 
as expected. The effective correlation function g^i^) can be expressed in 
terms of the width-correction A of the size distribution so that 

(1 + .4) -^(1^.4/8) 

^ — — ' ^ ^ 

with A = (a)^/(a^). Note that A is well defined for any size distribution 
function, so that Eq. ( pT[ ) can also be applied to polydisperse situations. In 
the limit of small volume fraction — > 0, one can estimate the normalized 
pressure by 

P,^(l + r)ug2aiiy)^^ , (18) 



12 Stefan Luding and Oliver Straufi 



as proposed by Zhang et al. [ pO[ , when disregarding the dependence of g{r) 
on the types of the coUision partners. The values of P/Pi in the hmit h' ^ 
agree very well with the simulations. Using the effective particle correlations, 
one can define 

P2iiy)^^-l = il+ r^gAi^) , (19) 

and compare the resulting expected reduced pressure with the simulation 
results from Fig. ||. An almost perfect agreement between P and P2{v) is 
obtained for v < 0.4 and even for larger ly « 0.65, the difference is always less 
than about two percent, and, the quality factors for all simulations collapse. 
Note that the quality is perfect (within less than 0.5 percent for all < 0.65) 
if P2{i^) is multiplied by the empirical function 1 — i^^/10, as fitted to the 
quality factor P/P2- Thus, based on our simulation results, we propose the 
corrected, nondimensional mixture pressure 

Piiiy) = ^-l = il + r),y9A{,^)[l-agiy^] , (20) 

with the empirical constant Ug « 0.1, for the pressure for all v < 0.65. 
For larger v the excluded volume effect becomes more and more important, 
leading to a divergence of P/ Pa- Furthermore, in the high density regime, the 
behavior is strongly dependent on the width of the size distribution function, 
see Fig. |^. 



3.4 Accounting for the dense, ordered phase 

The equation of state in the dense, ordered phase has been calculated by 
means of a free volume theory [^^,|3^, that leads in 2D to the reduced 
pressure Pfv = l/(\/tWx7i^~ 1) with the maximum volume fraction J^max- 
Based on our simulation results we propse the corrected high density pressure 

/'de„se = ^ [l+adiVra^.^-vy] , (21) 

where the term in brackets [. . . ] is a fit function with — 0.340 and Cp = 
1.09. The special case = leads to the theoretical result Pfv. In the left 
panel of Fig. ^, data from set B (see table ||) are presented, together with 
the "low" and "high" density predictions P4 and Pdcnsc, respectively (dashed 
and dotted lines). 

To our knowledge, no theory exists, which combines the disordered and 
the ordered regime. Therefore, we propose a global equation of state 

Q = P4 + m{v) [Pdcnsc - Pa] , (22) 



with an empirical merging function 



Poly disperse granular gases 13 



0.95 



0.85 



1.01 



8." 



0.99 



ax 

1 o o 



0.4 0.8 

J I L 



0.2 0.4 o.e 
V 



0.8 



Pi/Po 

Pi'Po 
Wq=0 

Wo=0.25 

Wo=0.50 

R=0.50 

R=0J5 



Fig. 7. Quality factor P/Po from the inset of Fig. ^. The lines give P2(v)/Po from 
Eg. In the inset, the simulation data for P axe rescaled by Piiv) from Eq. 



which selects P4 for v Vc and Pdcnsc for v ^ Vc with the width of the 
transition ttiq. In Fig. ^, the fit parameters Vc ~ 0.70 and mp ~ 0.009 lead to 
qualitative and quantitative agreement between Q (thick line) and the simu- 
lation results (symbols). However, also a simpler version Qq (thin line) with- 
out numerical corrections leads to reasonable agreement when mg = 0.015 is 
used, except for the transition region. The pressure drop when v is increased 
above v,. is qualitatively reproduced but no negative slope occurs. Due to the 
latter fact, the expression Qq allows for an easy numerical integration of P. 
We selected the parameters for Qo as a compromise between the quality of 
the fit on the one hand and the treatability of the function on the other hand. 

Remarkably, as one can see from Fig. |^ (Right), the dimensionless pres- 
sure Q from Eq. ( |2^ ) describes, at least qualitatively, the behavior of the 
polydisperse simulations when f,„ax = 0.858 is used. Note that the pressure 
drop at the transition Vc ~ 0.7 from the low density, disordered regime to the 
high density, ordered regime, is almost non-existent, since P4 and Pdense are 
almost collapsing in this range of density. 

4 Pressure gradient due to gravity 

In an experiment on earth, usually gravity plays an important role, it in- 
troduces a pressure gradient. Therefore, the density and pressure profiles of 
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Fig. 8. Dimensionless pressure P from set B simulations (symbols) plotted against 
the volume fraction v. (Left) The dashed lines are P4 from Eq. (^^ and Pdcnsc 
from Eq . (|2l| ) . The thick solid line is Q, the corrected global equation of state from 
Eq. ^) with the fit-parameters Ug = 0.1, = 0.340, ap = 1.09, = 0.701, 
i-'max = 0.9069, and mo = 0.00928. The thin solid line is Qo without corrections, 
i.e. Og = 0, Ud = 0, and mo = 0.0015 and Vc — 0.7, so that Qo = P2+m(i/)[Pfv— P2]. 
(Right) Analogous data from set Dl, where only J^max = 0.858 is different from the 
left panel. 



granular systems in equilibrium in the gravitational field are examined in the 
following. Here, a horizontal wall at z = is introduced in a periodic two di- 
mensional system of width L — lx/{2a), infinite height, and the gravitational 
acceleration g — 1 ms^^. 



4.1 Density profile in dilute systems 



In the special case of low density, one can use the equation of state of an 
ideal gas and express the pressure as a function of the energy density: 

p - I = nT , (24) 

with the number density n = N/V = n(z) — ^{z) / {tto^) and the "granular 
temperature" T = E/N va two dimensions. 

The gradient of pressure dp/dz compensates for the weight nmgLdz of 
the particles in a layer of height dz, so that 

dp dp dv 

— ^——^-nmg (25) 
dz dv dz 
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Separation of variables and the assumption of a constant temperature leads 
to the density profile for an ideal gas 

/ N ( "fngiz ~ zo)\ i^Q 

= 1^0 exp I — I or z[i^) = zq + zt in — , (26) 

with v < h'o and zt = T/{mg). In a system with a constant particle number 
A^, one has 



iV = / v{z)dz = — I (z(z/) - 2o) di^ ■ (27) 



Eq. ( p7| ) allows to determine analytically the volume fraction i/j^ at the bottom 
zq, in the dilute limit, by integration of z(z^) 

defined here for later use. 

This case can be extended to dilute and weakly dissipative systems, since 
the temperature is almost constant except for the bottom boundary layer 
[^^^. In the following we rather extend it to arbitrary density, but keep 
r = 1. 



4.2 Density profile for a monodisperse hard sphere gas 

In the dense case, Eq. ( |2^ ) is modified to 

p^nT[l + 2vg2a{v)] , (29) 
using Eq. (psl) with r = 1, and inserting Eq. (|29|) into Eq. (psl) leads to 



^^Tz-d-J^^^' + ^'^^^^^'^^^Tz--^'^'- ^''^ 



Assuming again that T is constant, one gets 



diy ^ 



T |1 + - (2. I = T 3(,„^)3 , (31) 



which, inserted in Eq. (|30|), allows integration from to v and from zq to z: 

' ' ' ' -5!!^ 



• n I"' 1-"' (1 -■'')' (1 -"■)'] J' J„ 

and leads to an implicit definition of v(z): 



= ^'t^. (33) 

ZT 
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We express z as a function of the volume fraction 

= In In h 2g2a[vo) ~ ^92a[v) , (34) 

zt v a 1 — V 

with the unknown volume fraction at zq, which, however, is determined 
using Eq. (|^: 

r- ^^d= dv = Z/o --; r2 , (35) 

ztL Jo zt 8(1 - v^Y 

where A'o is the number of particles above a given height zq. Only if zq = 0, 
one has iVo — N . This leads to a third order polynomial for i/q, 

4 - ^^dvl + (16i^d + 8)z^o - 8i/d = , (36) 

which can be solved analytically |3^ , and always has at least one real solution. 
Note that the function gia(y) is wrong at high densities z/ > z^c, so that also 
the pressure is not correct for high densities. This fact is discussed also by 
D. Hong who performed the three dimensional calculations analogous to 
our 2D calculus in this section. 



4.3 Comparison with simulations 

In this subsection, the theoretical density profile in Eq. (|3^), with the param- 



eter z^o determined via Eq. (36), is compared to numerical simulations with 
the parameters as specified in table In Fig. ||, the rescaled height z/zt 
is plotted against the volume fraction z^, according to Eq. (|3^). Note that 
even when the simulation parameters are rather arbitrary, the data follow a 
master-curve from z^ = to z^ = z^o (or equivalently from z = oo to z = 0) 
only shifted vertically such that z(i'o) = 0. The agreement between simula- 
tion and theory is almost perfect, except for simulation IV where densities 
above v k, 0.65 are observed, i.e. above the limit of validity of the equation 
of state. Therefore, the numerical values of zj zt are systematically smaller 
than the theoretical line obtained for vq from table |[ 





N 


L/(2a) 


T (kgm^s-^) 








I 


1562 


100 


3.07 X 10~" 


29.4 


0.418 


0.240 


II 


3000 


100 


2.22 X 10"** 


21.2 


1.110 


0.396 


III 


1000 


100 


2.61 X 10"^ 


2.49 


3.151 


0.567 


IV 


1000 


10 


6.13 X 10"'-' 


5.85 


13.41 


0.755 



Table 4. Simulation parameters for density profile measurements. In these simula- 
tions, the particle radius a = 5 x 10~* m and the particle mass m = 1.047 X 10"" kg 
were not changed. 
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Fig. 9. Rescaled height zJzt plotted against the volume fraction v from four differ- 
ent simulations, see table 0, with logarithmic horizontal axis. Symbols are simulation 
data and the lines correspond to Eq. (H) with the respective value of vq. 



The only way to get the correct theoretical density profile is for points 
with V < 0.65. The value of zq = 13.5zt, where v(zo/zt) = 0.65, is taken 
from the simulation data and the normalization accounts only for the A'o 
particles above zq. The resulting density profile (dashed line) nicely agrees 
with the simulations, and its limit of validity is inidcated by the angle at 
z/zt = 13.5 and v = 0.65. 

The reason for the increased density at low z/zt in the case of simulation 
IV is the pressure drop due to crystallization in the equation of state, see Fig. 
0. A higher density is necessary to sustain a given pressure when v > 0.65. 
The data for the lowest z/zt « are slightly off due to the wall induced 
ordering at z/zt — 0. 

If, instead of 2vg2a(v) in Eq. (30), we use the more general form Qo, we 
have to integrate the differential equation dp/dz = vmg j (tto?) numerically 
with p = nT{l + Qo) and the condition that Eq. ( p7| ) is fulfilled. Simula- 
tion IV and the numerical solution are compared in Fig. ^ The qualitative 
behavior of the density profile is well reproduced by the numerical solution 
with i/(zo = 0) = 0.8016. Note that the averaging result is dependent on 
the binning - we evidence strong coarse-graining effects in the dense, ordered 
region with densities v > 0.70. 
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Fig. 10. Rescaled height z/zt plotted against the volume fraction u from simulation 
IV, see table ^ The symbols are simulation data of the particle-center with a rough 
binning Az = 0.01007 m. The solid line is the density from the same data but 
with a much finer binning Az = 0.00184 m. Both binnings start at zo + tci with 
zo = 0. The dashed line corresponds to the numerical solution of Eq. (^) with 
p — nT{l + Qo)- The transition density z/c ~ 0.7 is indicated by an arrow. 

4.4 Bidisperse systems with gravitation 

In Fig. |ll| the species volume fractions i^i and 1^2 are plotted against the 
vertical coordinate z/(2ai). The data are obtained after long equilibration, 
in a system with N = 2000 particles, width ~ Ix/i'^ai) — 100, and the 
size distribution of set Dl in the previous section, with ai = 0.0005m. The 
particles are no longer mixed, as in the homogeneous, periodic systems of 
the previous section. Segregation takes place, the larger and heavier particles 
(squares) settle close to the bottom, whereas the gas of small and lighter 
particles (circles) extends to larger heights. Knowing both volume fractions, 
one could compute ^ as a function of the height, and insert it into Eq. ( pT|) 
in order to compute the density profile. 

5 Summary and Outlook 

In summary, we reviewed existing theories for dilute and dense almost elas- 
tic, smooth, 2D granular gases. For mono-, bi-, and polydisperse systems, 
we compared theoretical predictions with numerical simulations of various 
systems. The collision frequency, the energy dissipation and the equation 
of state, i.e. the scaled pressure, are nicely predicted by the theoretical ex- 
pressions up to intermediate densities. Especially, for arbitrary particle size 
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Fig. 11. Dimensionless height 2/(2ai) plotted against the volume fraction (in loga- 
rithmic scale). The thick, dashed line gives 1 — A (which is zero when only the small 
species is present at large heights), while the solid line indicates the total volume 
fraction of the mixture i^™ = 1^1 + ^2. 



distribution functions, the equation of state can be written in a nice form 
which only contains the width-correction A of the size-distribution function. 
A small, empirical correction can be added to the theories to raise the qual- 
ity even further. Finally, a merging function that connects the low and high 
density theories is proposed to give a global equation of state for all densities 
and size distribution functions. 

The equation of state is used to compute analytically and numerically the 
density profile of an elastic, monodisperse granular gas in the gravitational 
field. When a mixture is simulated, segregation is observed, a case to which 
the theory cannot be applied. For densities below v ~ 0.65, the analytical 
solution works well, for higher densities close to the maximum density, one 
has to use a numerical solver, since the global equation of state cannot be 
integrated analytically. The strange shape of the density profile, as obtained 
from simulations, is nicely reproduced. 

The simulations and the theories presented here were applied to homoge- 
neous systems. The range of applicability may be reduced by the fact that 
already weak dissipation can lead to strong inhomogeneities in density, tem- 
perature, and pressure. In a freely cooling system, for example, clustering 
leads to all densities between z/ « and v « fmax- The proposed global equa- 
tion of state is a necessary tool to account for such strong inhomogeneities 
with very high densities, above which the low-density theory fails. For another 
approach to handle the high density regions, see Ref. p5[ |. 
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The proposed global equation of state is based on a limited amount of 
data. It has to be checked, whether it still makes sense in the extreme cases of 
narrow w{a), where crystallization effects are rather strong, and for extremely 
broad, possibly algebraic w{a), where A is not defined. What also remains 
to be done is to find similar expressions not only for pressure and energy 
dissipation rate but also for viscosity and heat- conductivity and to extend 
the theory to three dimensions 
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